clear
load('./data/initial.mat');
ts = 70;
te = 80;
vx = zeros(nx,ny);
Nu = zeros(nx-1,te-ts+1);
for j = ts:te
    load([saldir 'dat',sprintf('%4.4d',j)]);
    vx(2:nx-1,2:ny-1)=-dy1p5*(phii(2:nx-1,3:ny)-phii(2:nx-1,1:ny-2));
    hfconv = mean(vx(:,2:end-1).*Ti(:,2:end-1),2);
    hfcond = -0.5*dif*diff(T0i)/dx;
    hf1 = 0.5*dif*(T0i(1,10)-T0i(end,10));
    Nu(:,j-ts+1) = (hfconv(1:end-1)+hfcond(:,10))./hf1;
end
Numean = mean(Nu,2);